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INTRODUCTION  AND  BACKGROUND 


One  of  the  most  interesting  phenomena  seen  in 
satellite  images  in  the  visual  wavelengths  is  sunglint, 
sometimes  referred  to  as  sun  glitter.  Over  oceanic  regions  or 
even  lakes  and  rivers  in  the  tropic  and  subtropic  regions,  there 
appear  to  be  areas  of  increased  brightness  or  reflectance  in  the 
observed  image  where  there  are  no  obvious  clouds.  This  is  the 
phenomenon  of  sunglint,  which  is  caused  by  the  reflectance  of  the 
sun  from  a  "shimmering"  water  surface  at  just  the  correct 
relative  sun-observer  position  relationships.  The  slopes  and 
directions  of  the  wave  facets  are  such  that  the  reflected  solar 
energy  is  directed  into  an  observing  system,  such  as  a  satellite 
sensor.  The  original  studies  of  sungl int/wind  relationships  were 
done  from  an  aircraft. 

Locating  sunglint  patterns  is  a  problem  composed  of 
two  parts.  The  first  is  to  locate  the  sun  and  satellite  on  a 
universal  coordinate  system,  as  well  as  to  calculate  all 
geometrical  quantities  and  relationships  as  a  function  of  date 
and  time.  This  is  basically  a  problem  in  spherical  trigonometry. 
The  second  part  is  a  determination  of  the  relationships  as  a 
function  of  various  angles  of  reflectance  and  the  sea  state 
conditions.  The  NEPRF  system  is  based  on  the  work  of  Cox  and 
Munk  (1954). 

In  1979,  NEPRF  acquired  a  computer  program  to  calculate  the 
sunglint  pattern.  This  program,  coded  to  run  on  the  Control 
Data  Corporation  (CDC)  6500,  produced  a  hard  copy  on  the  Varian 
printer.  Isopleths  of  reflectance  were  produced  for  a  constant 


wind  speed  in  the  sunglint  region.  The  final  graphic  chart  was 
placed  over  hard  copy  DMSP  images  using  a  light  table.  The 
image  and  graphic  were  co-registered  by  using  the  technique  of 
Tsui  and  Fett  (1980).  The  latitude  and  longitude  of  the  primary 
spectral  point  (PSP)  were  determined  by  hand,  and  the  scales  of 
the  graphics  and  the  image  were  made  the  same.  The  system  was 
completely  manual.  Input  into  the  reflectance  program  was  also 
done  manually,  with  the  required  parameters  such  as  satellite 
height,  inclination  angle,  time  and  longitude  of  the  ascending 
node,  and  the  time  of  the  image  entered  on  punched  cards. 

With  the  advent  of  relatively  inexpensive  mini/micro 
computers  and  graphics,  such  as  the  Hewlett-Packard  (HP)  9020 
and  the  Zenith-248,  it  is  now  possible  to  display  the  satellite 
images,  along  with  the  overlays  of  reflected  intensities  and  the 
navigation  grids  and  geography,  with  relative  ease.  Input  to  the 
new  system  is  automatic;  an  ephemeris  model  is  a  part  of  the 
input.  Parameters  may  be  obtained  from  regular  messages  or  other 
computer  files  available  to  the  user. 

The  original  NEPRF  sunglint  program  was  developed  from 
code  written  by  Irwin  Ruff  of  NOAA-NESS  during  the  late  1970's. 
His  code  was  modified  for  the  HP-9020  using  Fortran  77,  along 
with  the  graphics  processing  and  operating  system  developed  by 
Global  Imaging,  Inc.  NEPRF  never  received  any  program  documenta¬ 
tion,  physical  description  or  any  other  reference  to  Ruff's 
model.  The  in-line  comments  that  he  placed  in  the  code  were  very 
good,  however.  Using  these  in-line  comments  and  some  of  NEFRF's 
own  modifications,  this  report  will  describe  the  inner  workings 
of  the  Ruff  (NEPRF)  model. 
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2.  THE  RUFF  MODEL 

At  this  point  in  the  model  description,  we  will  assume  that 
all  required  sunglint  input  parameters  are  available.  The  input/ 
output  methodologies  will  be  described  in  Section  3. 

To  perform  field  calculations  of  any  geophysical  parameter, 
a  computational  grid  or  coordinate  system  must  be  devised.  Ruff 
chose  to  use  the  Satellite  Normal  Universal  Grid  System  (SNUGS). 
In  this  system  the  positive  X  axis  is  directed  along  the  sub¬ 
satellite  track;  the  Y  axis  is  directed  along  the  satellite  scan 
lines  and  is  positive  to  the  right  of  the  track.  Ruff  used  the 
names  "orbital  longitude"  and  "orbital  latitude,"  respectively, 
to  name  his  coordinates.  The  sun's  position  was  the  "solar" 
orbital  latitude  and  longitude.  The  origin  of  the  system  is  the 
equator  and  the  longitude  of  the  ascending  node.  The  region  for 
calculations  is  defined  by  the  sun's  position  in  the  along-track 
path  and  by  the  satellite's  across-the-track  scan.  The  grid  is 
illustrated  in  Figure  1.  An  increment  of  one  degree  is  used 
along  the  track,  while  the  increment  across  the  track  is  five 
degrees  of  satellite  scan  angle.  The  array  size  is  61  by  17. 

With  an  orbital  period  of  100  minutes,  the  time  increment 
represented  is  approximately  17  minutes. 

The  first  task  in  the  program  is  to  determine  the  approxi¬ 
mate  solar  hour  angle  (the  sun's  ground-point  longitude  at  the 
zenith),  the  relative  angle  between  the  sun  and  satellite,  and 
the  solar  declination.  Given  the  date  and  the  time  of  the 
satellite's  ascending  node,  the  declination  IS)  is  calculated 
as  follows.  A  sine  curve  is  fitted  to  the  time  differential, 
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Figure  1.  The  computational  grid, 
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Del(DTG),  between  OOZ  on  March  21  and  the  date  time  in  hours 
(plus  tractions)  of  the  satellite's  ascending  node,  using  an 
amplitude  of  23.5  degrees.  The  argument  (ARG)  of  the  sine  curve 
is  as  follows: 

ARG  *  Del(DTG)/24. 0/365. 0*360. 0* (PI/180 . 0) ,  and 
<f  «  23 . 5*sine ( ARG) 

The  solar  hour  angle  (HA)  is  calculated  as  follows: 

HA  =  540.0  -  Time (hours ) *15 . 0, 

where  Time  =  Hours  +  Minutes/60.0  +  Seconds/3600.0  and  HA  is 
limited  to  the  range  between  0  and  360  degrees.  Remember  that 
the  hour  angle  is  calculated  from  nadir;  that  is,  at  1200Z,  the 
sun  is  overhead  at  zero  degrees  longitude. 

The  solar  hour  angle  is  next  converted  to  a  local  hour  angle 
(LHA) .  The  LHA  is  the  angular  difference  between  the  longitude 
of  the  ascending  node  and  the  hour  angle  of  the  sun,  measured 
eastward.  All  calculations,  when  referring  to  longitude,  refer 
to  this  relative  angle,  the  LHA. 

Having  the  position  of  the  sun  in  geocentric  coordinates 
(declination  and  LHA),  we  now  transform  them  to  SMUGS  coordinates 
to  simplify  all  further  calculations.  Figure  2  illustrates  the 
geometrical  relationships. 

Given  the  declination  IS),  the  LHA  (X),  and  the  inclination 
of  the  satellites  equatorial  crossing  (i),  the  solar  orbital 
latitude  (4)  and  the  solar  orbital  longitude  (t)  may  be 
calcul?.ted  in  the  following  manner.  Following  the  geometry  of 
Figure  2,  given  spherical  triangles  (a)  EPS,  ESO,  ENS,  ESA,  use 
the  law  of  sines  for  a  spherical  triangle: 
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Slnx'  Sin  w  Sin(e)  +  l • )  sin  90° 

For  a  EPS.  -  =  -  ;  For  a  EBA,  - - -  - 

Sin  Q  Sin  <T  Sin  X'  Sin  Q 

Sin(u>+i '  )  Sin  90°  Sin  e  Sin  t 

For  A  ENS,  -  =  -  ;  For  A  ESO,  -  =  -  . 

Sin  *  Sin  Q  Sin  *  Sin  Q 

Angles  PAE,  ENO,  OEN,  PEA  are  right  angles.  Further, 
a)  +  e  +  i’  =  90°;  i*  =  90°  -  i;  **  =  90°  -  ♦  ;  Cos  90°  =  0.0; 

e  =  90°-(o)+i');  Sin  (e+i')  =  Coscu;  <J‘  =  90°-<f;  Sin  90°  =  1.0; 
e  +  1 '  =  90°  -  <u;  Sine  =  Cos (u-f i '  ) ;  Sin<P  =  Cos<J;  X'  =  X  -  360°; 

Cos  1'  =  Sin  1  ;  Sin*'  =  Cos*;  X*  =  -X  . 

For  a  ESA,  Cos  *  Sinr  =  Cos(o)  +  i')  Sin  Q.  Substituting 

for  Sin  Q  in  a's  EPS  and  EBA, 

Sin  id  Cos  <d  Cos  6  Sin  X' 

-  _  -  #  or  Tan  id  *  -  . 

Sintf'  SinX'  Sin  <f  Sin  <f 

Substituting  for  e  in  A  ESO  and  expanding, 

Cos  *  Sin  r  =  (Cos  a)  Cos  i'  -  Sin  id  Sin  i’)Sin  Q. 

After  substituting  for  Sin  (e  +  i * )  in  A  EBA,  now  substitute 

for  Sin  Q  in  the  expanded  form  of  a  ESO: 

Cos  *  Sin  r  =  (Cos  u>  Cos  i'  -  Sin  id  Sin  i')Sin<f/Cos  u>,  or 

Cos  *  Sin  r  =  (Cos  i'  -  Tan  id  Sin  i'JSlntf. 

Substituting  for  Tan  w.  Sin  i'  and  Cos  1',  we  obtain 

Sin  t  =  (Sin  1  Sin  6  +  Cos  i  Cos  i  Sin  X)/Cos  *. 

For  the  solar  orbital  latitude,  expand  a  ENS: 

Sin  *  =  (Sin  <d  Cos  i'  +  Cos  u>  Sin  i')  Sin  Q. 

Now  substitute  for  Sin  Q  from  a  EBA,  and  divide  by  Cos  u>;  we  then 

obtain. 

Sin  *  *  (Tan  <d  Cos  i'  +  Sin  i’JSin  6. 
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Now  substitute  for  Tan  u>,  Cos  i’,  and  Sin  i',  producing 

Sin  *  =  Cos  i  Sin  <f  -  Sin  i  Cos  <f  Sin  X. 

The  next  task  of  part  1  (see  Introduction)  in  the  initial¬ 
ization  process  is  to  calculate  the  geodetic  latitude  and 
longitude  of  the  starting  point  of  the  grid.  The  starting  point 
is  30.0  degrees  of  solar  orbital  longitude  prior  to  the  solar 
position  at  the  time  of  the  satellite's  ascending  node.  This 
geodetic  position  is  determined  by  an  integration  backwards  from 
the  equator  and  the  longitude  of  the  ascending  node,  with  an 
azimuth  equal  to  the  satellite's  inclination  angle  (i)  +  90.0 
degrees.  The  number  of  increments  used  in  the  integration  is 
30.0  minus  the  solar  orbital  longitude  of  the  sun.  For  example, 
if  the  solar  orbital  longitude  is  4.3,  the  integration  is  done  in 
25  one-degree  steps,  plus  one  step  of  .7  degrees.  The  purpose  of 
this  step  is  to  be  able  to  relate  the  calculations  on  the  grid 
back  to  the  image.  This  step  is  not  necessary  to  the  reflectance 
determination  and'  was  not  part  of  the  Ruff  model.  The  geodetic 
latitude  and  longitude  values  are  determined  at  the  same  time  the 
SNUGS  coordinates  are  calculated.  A  simple  method  is  used,  as  a 
high  degree  of  accuracy  is  not  required.  We  are  assuming  that  we 
are  looking  at  a  "snapshot''  with  only  a  limited  time  involved. 

The  integration  of  position  proceeds  by  providing  a  known 
latitude,  longitude  and  azimuth,  and  then  calculating  a  new 
latitude,  longitude  and  azimuth,  assuming  a  stationary  earth. 
Figure  3  Illustrates  the  geometry  of  the  integration. 
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Figure  3 
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The  Integration  o£  the  satellite's  position. 
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Using  the  law  of  cosines  and  sines  for  spherical  triangles, 
Cos  *'  =  Cos  St  Cos  +  Sin  tfr  Sin  * '  Cos  it,  where 

*'  =  90.0  -  *;  Cos  *'  =  Sin  Q;  Sin  *'  =  Cos  *. 
Substituting,  we  obtain 

Sin  *  =  Cos  St  Sin  +  sln  Cos  *t  Cos  *t*  Similarly, 

Cos  =  Cos  *t  Cos  *t+l  =  Cos  “  Sin  *t  s*n  *t+-l/ 

or  rearranging, 

Cos  =  (Cos  St  -  Sin  *t  Sin  *t+i)/(Cos*t  Cos  *t+l)* 

And, 

Sin  it  +  l  =  Cos  +  ^  Sin  i^/Cos 

The  program  is  now  ready  to  make  reflectance  calculations  on 
the  grid,  but  additional  geometrical  relationships  are  needed 
first.  As  described  previously,  the  scan  line  increment  is  five 
degrees,  and  the  along-the-track  increment  is  one  degree  of  solar 
orbital  longitude.  We  start  40.0  degrees  to  the  left  of  the 
track  and  30.0  degrees  from  the  solar  orbital  longitude  of  the 
sun.  The  reflectances  are  first  determined  across  the  track, 
then  along  it.  The  first  quantities  needed  at  the  different  scan 
angles  are  the  zenith  angle  of  the  satellite  and  its  geocentric 
angle . 

Figure  4  illustrates  the  geometry  between  the  various 
angular  measurements.  The  law  of  sines  for  a  plane  is  used  to 
relate  them: 

Sin  ©'  *  (r  +  h)Sin  h/r,  and 
e1  =  180.0  -  ©,  c  ■  B  -  h,  where 

9  *  the  zenith  angle  of  the  satellite, 

h  =  the  scan  angle  of  the  satellite. 
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<7  ■  the  geocentric  angle, 

h  =  the  height  o£  the  satellite,  and 

r  -  the  radius  of  the  earth. 

The  zenith  angle  of  the  sun  (?)  is  shown  for  reference. 

With  the  satellite  geocentric  angle  (<r),  along  with  the 

solar  orbital  latitude  (*)  and  longitude  (t),  the  zenith  angle 

(?)  and  its  azimuth  (a)  to  the  sun  at  the  particular  grid  point 
are  acquired.  The  geometry  for  the  acquisition  is  shown  in 
Figure  5. 

From  the  spherical  triangle  VSO,  the  zenith  angle  can  be 
written  directly  using  the  law  of  cosines: 

Cos  ?  *  Sin  *  Sin  <7  +  Cos  *  Cos  O'  Cos  r. 

The  azimuth  requires  a  little  more  work. 

As  before,  o  =  90.0  -o';  a'  =  180.0  -  K;  *  *  90.0°  - 

For  a  VSE,  Cos  Q  =  Cos  ?  Cos  o  +  Sin  ?  Sin  o  Cos  a'. 

For  A  ESO,  Cos  Q  *  Sin  *  Cos  90°  +  Cos  *  Sin  90°  Cos  r  or, 

Cos  Q  =  Cos  *  Cos  T. 

Substituting  for  Cos  Q  in  h  VSE, 

Cos  ?  Cos  o  +  Sin  ?  Sin  o  Cos  a'  *  Cos  ?  Cos  o/Sin  o. 
Rearranging  and  dividing  by  Sin  o, 

Cos  o'  Sin?  =  Cos  *  Cos  r/Sin  <7  -  Cos  ?  Cos  <7/Sin  <7. 


Substituting  for  Cos  ?, 

2 

Cos  ♦  Cos  t  Cos  *  Cos  <7  Cos  r 

Cos  «'  Sin  ?  =  -  -  Sin  ♦  Cos  <7  -  - 

Sin  v  Sin  <7 

or 

Cos  *  Cos  r  _ 

Cos  ot'  sin  ?  *  -  (1.0  -  Cos  «)  -  Sin  *  Cos  <7, 

Sin  v 


or 

Cos  a  *  Sin  *  Cos  <7  -  Sin  <7  Cos  *  Cos  t. 
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Figure  5.  The  relationship  of  azlsuth  and  polar  zenith  to  SMUGS 
coordinates . 
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The  geodetic  coordinates  are  now  generated.  With  these 
latitudes  and  longitudes  the  surface  wind  field  can  be  extracted 
from  the  FNOC  global  marine  data  fields  for  use  as  a  variable 
wind  speed  in  the  reflectance  calculations;  the  graphics  and  the 
image  can  then  be  co-registered. 

All  the  parameters  necessary  for  the  determination  of  the 
solar  reflectance  are  now  present. 

There  are  three  steps  to  determine  the  intensity  of  the 
reflected  solar  radiation.  First,  the  reflectance  coefficient 
(P)  which  is  only  a  function  of  the  incidence  angle  (reflected 
angle),  and  is  an  approximation  of  the  Fresnel  reflection 
coefficient,  must  be  computed.  The  form  of  the  Fresnel 
reflection  coefficient  is: 

P  =  Cos  P  ( (exp(-kl*Cos  P))  +  k2*Cos  P), 
where  kl  *  6.09  and  k2  =  .019. 

The  angle  of  reflection  (the  incidence  angle),  beta  is  the 
angle  necessary  to  reflect  incoming  solar  radiation  in  the 
direction  theta,  with  an  azimuth  angle  alpha.  Figures  6a  and  6b 
illustrate  the  geometry  of  the  problem.  For  spherical  o  TVS, 

Cos  2P  *  Cos  9  Cos  ?  +  Sin  9  Sin  ?  Cos  a,  or 

P  =  Cos-1  ( ( 0 . 5* ( 1 . 0  +  Cos  2P))1/2). 

The  second  step  is  determining  the  probability  that  the 
slopes  and  directions  of  the  wave  facets  of  the  sea  surface  will 
reflect  the  Incoming  solar  radiation  at  the  correct  angles  to 
give  perfect  reflection.  This  probability  (P)  is  a  function  of 
the  variance  of  the  wind  speed  and  the  slopes  of  the  sea  surface: 

P  *  (exp (-Tan2  e?/<r2)  -  exp(-Tan29j!/v2 ) )  (a?  -  aj )/2rr,  where 
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Figure  6.  The  relationship  between  satellite  zenith  angle/ 
solar  zenith  angle/  aziauth  and  Incidence  angle. 

(a)  Plane  view,  (b)  Section  view. 
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subscript  i  and  M  refer  to  lower  and  upper  limits,  respectively, 
of  the  Integration.  Superscript  N  refers  to  the  fact  that  we  now 
are  measuring  angles  relative  t-o  a  normal  to  the  sea  surface. 


This  probability  is  an  integration  over  a  finite  solid  angle 
<fw.  The  limits  are  between  ©i,  aj  and  ©£,  o$  centered  on 
©c •  ac*  Tan  ®c  i®  the  slope  of  the  sea  surface  (wave  facet) 
or  the  tangent  of  the  zenith  angle  of  the  surface  normal  required 
to  reflect  at  a  given  6  and  azimuth  a.  Tan  ©£  is  the  ratio  of 
the  sine  of  the  vector  difference  of  e,  :  and  the  cosine  of  the 
vector  sum  of  ©,  <?: 


Tan  = 


Sin  ©c 
Cos  ©c 


Using  the  law  of  cosines  for  a  plane,  the  sine  and  cosine  of  the 
slope  normal  are  calculated  in  the  following  manner: 

Sin2  ©”  =  Sin2©  +  Sin  2 f  +  2*Sin©  Sin  ?  Cos  a”'  and 


Cos  ©c  *  Cos  ©  +  Cos  <T 


.N 


The  corresponding  azimuth,  ac  is  calculated  as  follows.  Again, 
using  the  law  of  cosines  for  a  plane. 

Sin2©  =  Sin2©”  +  Sin2-:  -  2*Sln©”  Sin  :  Cos  oj'  . 
Substituting  for  Sin2©, 

Sin2©”  *  Sin2©”  +  Sin2-:  -  2*Sin©”  Sin  i  Cos 

+  Sin2-:  +  2*Sin  ©  Sin?  Cos  a”,  or 
2*Sin2:  -  2*Sin©c  SinC  Cosa”  +  2*Sln©  Sin?  Cos  a”  *  0.0. 
Simplifying  and  rearranging, 

Cos  a”  «  (sin  :  +  Sin  ©  Cos  <x)/sin  ©”  • 

Figure  7a  and  7b  show  the  geometrical  relationships. 


16 


Figure 

(a 


The  vave  facets  slope  and  aziauth  relationships. 
Plane  view,  (b)  Section  viev. 
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For  the  solid  angle  integration  and  the  probability,  the 
limits  are  .01  radians  on  either  side  of  the  center  angles  ©jj, 

Kq  which  are  the  exact  slope  and  azimuth  needed  to  reflect  the 
ray.  Thus, 

ejl  *  ©c  +  “c  *  “c  +  <*“>  and 

©!-©?-  SS.  “t  *  a£  -  *<x. 

The  wind  dependence  is  introduced  into  the  probability  as 
the  variance  of  a  circular  Gaussian  distribution  of  the  normals 
of  the  slopes  of  the  sea  surface  (wave  facets).  This  dependence 

is  linear,  and  it  is  after  Cox  and  Munk  (1954): 

2 

o  =  A*v  +  B,  where 

2 

A  =  .0051,  B  =  .003,  w  is  the  wind  speed  in  m/sec,  and  <r  is 
the  variance. 

The  third  step  for  reflectance  is  to  determine  the  elemental 
solid  angle  over  the  same  limits  that  were  used  in  the  probabil¬ 
ity  statement.  The  integration  is  in  steps  of  .001  radians  for 
both  ©i  and  oc^.  The  increment  of  solid  angle  is  integrated 
on  a  Lambert  azimuthal  equal  area  projection  of  the  upper 
hemisphere,  centered  on  the  solar  direction  ©  =  i';  a  =  0. 

Figures  8a,  8b,  and  8c  show  the  relationships  between 
incidence  angle  and  azimuth,  as  referred  to  the  sphere  and 
the  projection.  Using  the  law  of  cosines  for  a  sphere,  we  may 
determine  the  incidence  angle  P}  for  &  TVS  in  the  solid  angle 
integration  (Fig.  8a): 

CosPj  *  Cos  Cos  ©^  +  Sin  ?  Sin  ©^  Cos  ajJ  . 
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Figure  8.  The  geoaetrical  relationships  for  solid  angle  calcula¬ 
tions  between  the  sphere  and  the  projection.  (a)  Plane  view 
of  sphere,  (b)  Section  view,  (c)  The  projection. 
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The  azimuth  from  the  sun  projected  point  ( M )  may  be  determined 
as  follows,  again  using  the  lav  of  cosines: 

Cos  ©i  *  Cos  Mi  Sin  Z  Sin  Pi  +  Cos  Pi  Cos  Z  . 

Substituting  for  Cos  Pi,  we  obtain 
Cos  ©j[  =  Cos  Mi  Sin  Z  Sin  Pi 

+  Cos  Z  (Cos  Z  Cos  eMi  +  Sin  Z  Sin  eNi  cos  aNi). 
Simplifying  and  dividing  by  Sin  Z,  we  produce 


Cos  €>i 

Cos  Mi  Sin  Pi  =  - 

Sin  Z 


Cos  Z  [cos  ai  Sin  ©i 


+ 


Cos  ©i 


Cos  Z ^ 
Sin  Z 


Cos  Mi  Sin  Pi  =  Cos  Z  Sin  ©Ni  Cos  «i  - 


Cos  ©'i 


Sin  Z 


(cos2-:  -  i.o). 


or 


Cos  Mi  *  (Cos  e?  Sin  Z  -  Cos  «:  Sin  ©*  Cos  oti)/sin  Pi 


The  solid  angle  is  calculated  on  the  projection  by 

determining  the  area  for  each  sector  and  combining  them  in  the 

correct  sequence.  <fu>  is  the  area  of  sectors  ABO  +  ADO  minus 

the  area  of  sectors  DCO  +  BCO.  The  area  of  each  sector  is  an 

integration  over  P  and  M  where 

ABO  V 
ADO  <" 

Areas  DCO  =  *  . 5*(2P)$Mi  / 

BCO  *“ 
x 

with  the  index  limits,  x  and  y,  a  function  of  each  sector. 
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The  Integration  of  solid  angle  begins  with  the  smallest 

values  of  and  0^*.  is  increased  first,  followed  by  0N£. 

N  N 

a i  is  now  decreased,  followed  by  0^.  In  this  way  the  proper 

signs  are  maintained  during  the  calculations. 

Throughout  the  program  there  are  numerous  instances  where 

anomalous  conditions  exist.  For  example,  the  zenith  angle  0  may 

|| 

be  equal  to  zero,  or  the  azimuth  angle  a  may  be  equal  to  zero. 

Ir  addition,  computations  must  be  made  using  double  precision. 

Once  the  integration  of  the  solid  angle  increment  is 
complete,  the  intensity  of  the  reflected  energy  for  any  point  on 
the  computational  grid  is  calculated  as  the  product  of  the 
•probability  (P)  and  the  reflectance  coefficient  (P),  normalized 
per  unit  solid  angle  (<foj).  Figure  9  shows  the  reflection 
coefficient  as  a  function  of  incidence  angle,  and  Figure  10 
depicts  the  probability  as  a  function  of  wave  slope  for  different 
wind  speeds.  The  reflectance  coefficient  and  probability  are  the 
two  terms  in  the  reflectance  (R)  calculations,  along  with  the 
normalization  factor  (the  solid  angle,  delta  omega): 

R  =  P  P/<fu). 

The  reflectance  is  converted  to  percent,  and  along  with  the 
grid  point's  latitude  and  longitude,  is  written  as  an  ASCII 
record.  Each  point  is  considered  as  a  randomly  spaced  data  point 
as  required  in  the  Global  Imaging  processing  system.  All  data 
with  reflectance  less  than  0.1\  are  set  to  zero.  After  the 
reflectances  have  been  determined  over  the  computational  grid. 
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Figure  10.  Probability  as  a  function  of  wave  slope  (degrees) 
and  wind  speed  (a/s). 
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the  system  will  perform  interpolations  to  place  the  reflectances 
to  a  regular  grid,  corresponding  to  the  same  lines  and  elements 
as  the  requested  image. 

3.  THE  NEPRF  SYSTEM 

The  capability  to  examine  images  with  the  reflectance 
patterns  superimposed  is  now  part  of  the  procedures  of  the  HP 
9020  graphics  and  image  processing  system.  The  only  input  to  the 
system  is  a  filename.  In  this  file  are  all  the  parameters  needed 
to  display  the  isopleths  of  reflectance  over  this  image  and  they 
are  automatically  obtained.  A  message  containing  the  orbital 
elements  is  read  and  decoded.  This  message  gives  the  epoch  date, 
time  and  the  inclination  angle,  the  argument  of  perigee,  the  mean 
anomaly,  the  eccentricity,  the  right-ascension  of  the  ascending 
node  and  the  mean  motion.  The  mean  motion  is  converted  to  the 
semi-major  axis.  These  parameters  are  fed  to  an  ephemeris 
program  that  gives  the  latitude,  longitude  and  height  of  the 
satellite  as  a  function  of  time.  The  time  interval  is  set  to  one 
minute.  The  ephemeris  file  is  scanned,  looking  for  the  ascending 
nodes.  The  closest  one  to  the  time  of  the  image  is  chosen.  The 
longitude,  time  and  height  of  the  satellite  are  extracted  on 
either  side  of  the  equator  and  then  are  linearly  interpolated  to 
the  equatorial  crossing  of  the  satellite.  There  is  one  more 
parameter  needed  to  run  the  sunglint  code.  This  is  the  wind 
speed.  There  is  an  option  to  give  the  program  one  constant  value 
of  wind  speed  (as  the  old  system  did),  or  to  extract  the  wind 
components  u  and  v  from  the  FNOC  marine  wind  field  analysis  at 
the  time  closest  to  the  image. 
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The  system  now  has  all  parameters  needed  to  make  the  sunglint 
reflectance  calculations,  navigate  them  and  display  the  results. 
This  information  is  passed  to  the  program  via  a  data  file.  After 
the  sunglint  module  has  been  run,  the  desired  image  is  displayed 
on  the  CRT,  followed  by  an  overlay  of  the  contours  of  reflectance. 
The  image/overlay  display  package  is  part  of  the  Global  Imaging, 
Inc.  processing  and  display  package  written  for  the  HP  9020. 

Their  system  manages  all  the  files  that  have  been  transferred 
between  the  SPCU,  the  Zenith  248  and  the  HP  9020.  At  this  time 
only  data  from  NOAA  AVHRR  imagery  can  be  processed  on  the  Global 
Imaging  System.  However,  DMSP  images  will  soon  be  added  to  the 
system's  capability.  A  difficulty  arises  from  the  navigation 
of  the  DMSP  satellites.  Another  method  of  running  the  system, 
especially  for  DMSP  images  from  the  Fleet  Numerical  Oceanography 
Center  (FNOC) ,  is  to  perform  the  sunglint  calculations  for  one 
12-hour  period  (seven  orbits).  In  this  way,  data  from  FNOC's 
specialized  data  base  could  be  extracted  and  navigated  for  very 
large  sectors  of  the  globe.  The  corresponding  reflectance 
overlays  would  be  present  for  display.  Figure  11  is  an  example 
of  the  display  of  the  reflectance  calculation  for  constant  wind 
speeds  (calm  and  5  m/sec) . 

At  this  time,  we  have  a  system  which  will  automatically 
display  an  image  along  with  the  corresponding  sunglint  patterns . 
The  wind  field  may  be  variable  or  constant,  at  the  user's 
discretion.  In  areas  of  the  world  where  the  sunglint  is  an 
observable  phenomena,  meteorologists  under  certain  conditions 


26 


Figure  11.  An  example  of  the  sunglint  pattern  for  a  constant 
field  of  wind  speed  (0.0  and  5.0  m/sec). 
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can  determine  in  a  qualitative  sense  the  correctness  of  an 
analyzed  wind  field.  Further  work  is  needed  to  determine  if 
a  quantitative  value  of  the  wind  may  be  determined  from 
observations  of  sunglint. 

Appendix  A  is  a  listing  of  the  sunglint  program,  with 
associated  subroutines;  Appendix  B  details  the  running 
instructions  for  the  sunglint  program.  Both  appendices  have 
been  developed  as  HP  9020  using  Global  Imaging,  Inc.  support 
software. 
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APPENDIX  A 

SONGLINT  PROGRAM  LISTING 


CCCNoption  trace  on 
PROGRAM  SNGLNT 
C  PROGRAMMER  IRWIN  RUFF 
C  PROGRAMMER  T.L.  TSUI 
C 

IMPLICIT  REAL  *8(A-H),  REAL  *8(0-Z) 

C 

COMMON  /ABC/  IDATE (3) , IPLOT, NUM1 , NUM2, NUM3, NUM4, NUM5 

REAL*4  PLAT, PLON, PLONT, PLATMAX, PLONW, PLATMIN, PLONE 

REAL*4  W(73, 144,2) , WIND 

INTEGER  TN, IQNX  (3) , LH 

REAL  *8  LN 

REAL  *4  DUMMY 

CHARACTER  CDATE* 12, CIQNX*12, CTN*4, CLH*12 
EQUIVALENCE  ( IDATE, CDATE) ,  < IQNX, CIQNX) ,  (TN, CTN) 

PI=DAC0S<-1.0D0) 

C 

C  INITIALIZE  GAE  PARAMETER  BLOCK  "sunalint" 

C 

CALL  initbl k ( "sungl int " ) 

C 

C*********************************************************************** 
C  INPUT  DATA  FROM  GAE  PARAMETER  BLOCK 
call  intgae (6h/YEAR/ , IYR) 
call  intgae(7h /MONTH/, IMO) 
call  intgae(5h/DAY/, IDAY) 
call  intgae (6h/H0UR/ ,  IHR) 
call  intgae(5h/MIN/, IMIN) 
call  real gae (8h/EQL0NG/ , DUMMY) 

LN=DBLE( DUMMY) 
call  strgae <7h/EQHEM/ , CLH) 
read ( clh, ’ (al ) ’ )  lh 
write<*,  '  <"EQHEM=  »',a2,)»)  LH 
call  realgaeOh/SATINCL/,  DUMMY) 

SIA=DBLE (DUMMY) 
call  real gae (Bh/SATALT / , DUMMY ) 

HT=DBLE( DUMMY) 

call  f ileop(6h/WIND/ , 5)  !  open  input  file  on  unit  5 

REWIND  5 
NW*1 
IPLOT =2 

write( cdate, *  (3i2.2) 9 )  iyr , imo, iday 
wr ite (CTN, *  (2i2. 2) 9 )  ihr,imin 
C 

20  READ (5, 5000, END *3999)  iend,jend 

READ(5, 5010, END *9999)  (((W(I,J,k),J=l, jend) ,1*1, iend) , k=l , 2) 

C 


c 

COLUMN 

FMT 

DESCRIPTION 

c 

IDATE 

1 

-  6 

A6 

YYMMDD  (DATE -TIME -GROUP) 

c 

c 

SI  A 

7 

-  10 

F4.  1 

SATELLITE  INCLINATION  ANGLE 
(E.G.  DMSP  F2:  SIA-98.7) 

c 

c 

HT 

11 

-  14 

F4.0 

AVERAGE  SATELLITE  HEIGHT  IN 
(E.G.  DMSP  F2:  HT-833) 

KM 

c 

c 

TN 

18 

-  19 

A4 

TIME  OF  LAST  EQUATOR  CROSSING  IN  ZULU 
(TN=HHMM,  HOURS  AND  MINUTES  OF  ZULU  TIME) 

20 


C  LN  20  -  24  F5. 1  LONGITUDE  OF  EQUATOR  CROSSING 

C  LH  25  A1  LN  LOCATION  INDICATOR,  E  -  EASTERN  HEMI . 

C  W  -  WESTERN  HEMI. 

C  IPLOT  27  II  >1  CONTOURS  AND  RADIANCE  VALUES  PLOTTED 

C  *1  CONTOURS  PLOTTED  ONLY 

C  <1  RADIANCE  VALUES  PLOTTED  ONLY 

C  NW  30  II  NUMBER  OF  CHARTS  DESIRED,  MAX.  6 

C  W  31-60  6F5.1  W(I)  IS  THE  WIND  SPEED  (M/SEC)  ASSOCIATED 

C  WITH  EACH  CHART. 

C*********************************************************************** 

c 


RADIUS  OF  EARTH 
R=6378. 16D0 

ORBITAL  INCLINATION  (EQUAL  TO  OR  LESS  THAN  90  DEGREES) 
ORBINC  *  180. DO -SI A 

ETA  IS  NADIR,  DELETA  IS  INTERVAL  IT  IS  READ 
DELETA=  5.0D0 

STARTING  POINT (LATITUDE)  ON  MAP  ENTERED  MANUALLY 
NCNT -0 

STARTING  POINT 
TAUI =0. ODO 

INTERVAL  REQUIRED 
DELTAU=1 . ODO 


C 

C 

C 

C  FIND  THE  SOLAR  DECLINATON 
C 

READ  (CDATE, 7001 )  IQNXC1) 

WRITE  (CIQNX, 7002)  IQNX(l) 

WRITE  (CDATE, 7003)  CDATE, TN 
IQNX ( 1 ) =IDIFDTG (CIQNX , CDATE) 

SO  L  =  FL  OAT ( IQNX ( 1 ) ) /24. DO/ 365 . D0*360 . DO*P I / 1 80 . DO 
S0LDEC=23. 5D0*DSIN (SQL ) 

C 

C  FIND  THE  LONGITUDE  OF  EQUATOR  CROSSING 
EQCX  =LN 

IFCLH.EQ. 1HW)  EQCX=360. DO-EQCX 


C 

C  FIND  THE  SOLAR  LONGITUDE 
C 

READ  (CTN, 7004)  SOL , AMDA 

S0L*S0L+AMDA/60. DO 

SOL  *DMOD ( 540 . DO -SOL  * 1 5 . DO , 360 . DO ) 

AMDA*EQCX -S0L+360.  DO 

AMDA  *DMOD ( AMDA , 360 . DO ) 

C 

WRITE  (6,6000)  IDATE 
WRITE  (6,6001)  HT, SI A 
IQNX ( 1 ) *1HE 

IF(SOL.LT. 180.)  GO  TO  12 
SOL *360. DO -SOL 
IQNX ( 1 ) *1HW 
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12  WRITE  (6,6002)  SOLDEC, SOL , IQNX ( 1 ) 

WRITE  (6,6003)  TN, LN, LH 
WRITE  (6,6004)  AMDA 
C 

c 

SNINC=SIN(ORBINC*PI / 180. DO) 

CSINC=C0S(QRBINC*PI/180. DO) 

SNDEC  =S  I N  ( SOLDEC*P I / 1 80 . DO ) 

CSDEC *CQS ( SOLDEC*P I / 1 80 . DO ) 

SNLMDA=SIN(AMDA*PI/1B0. DO) 

CSLMDA  *COS ( AMDA*P I / 1 SO . DO ) 

CCC  WRITE  (6,6006)  SN I NC, CS INC, SNDEC, CSDEC, SNLMDA, CSLMDA 

CALL  SLQRBC (SNINC, CSINC, SNDEC, CSDEC, SNLMDA, CSLMDA, SNSLOT, CSSLOT, 
1SL0BLG) 

CCC  WRITE  (6,6006)  SNSLOT, CSSLOT, SLOBLG 

C 
C 
C 

DO  11  K=  1 , NW 
C 
C 

WRITE  (6,6005)  K,WIND 
NUM1 =0 
NUM2  =0 
NUM3 -0 
NLJM4=0 
NUM5  =0 
C 

C  INITIALIZE  PARAMETERS  FOR  SUBROUTINE  SUBLTLN  -  -  SUBTRACK  LAT  LO 

NG 

C 

PHI  =  O.ODO 

AZ  =  SI A  -  90. ODO 

ALN  =  EQCX 

SNDTAU  =  DSIN(DELTAU*(dacos(-1.0d0)/180.0d0)) 

CSDTAU  =  DCOSCDELTAU* (dacos ( -1 . OdO) /ISO. OdO) ) 

C 

IF (NCNT. EQ. 0)THEN 
C 

TAUI =SL0BLG-30. ODO 

IFCTAUI. LT. O.ODO)  AZ *AZ+180. ODO 

BCKFRCT  =TAUI -INT (TAUI ) 

C 

C  STEP  A  FRACTION  OF  AN  INCREMENT  ALONG  SATELLITE  TRACK 

C 

IF ( ABS(BCKFRCT)  .GT.  0.0)  THEN 

SNDTAU  *  DS IN (BCKFRCT* (da cos ( -1 . OdO) /ISO. OdO) ) 

CSDTAU  =  DCOS (BCKFRCT* (da cos ( -1 . OdO) / 180. OdO) ) 

CALL  SUBLTLN (SNDTAU, CSDTAU, PHI , ALN, AZ , PHI2, ALN2, AZ2) 

SNDTAU  *  DSIN(DELTAU*(dacos ( -1 . OdO) / 180. OdO) ) 

CSDTAU  *  DCOS(DELTAU*(dacos( -1 . OdO) / 180. OdO) ) 

END  IF 

DO  100  1*0, INT (ABS (TAUI)) 

CALL  SUBLTLN (SNDTAU, CSDTAU, PHI, ALN, AZ , PHI2, ALN2, AZ2) 

PHI *PHI2 
ALN=ALN2 
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AZ  =AZ2 
100  CONTINUE 
END  IF 

IFCTAUI. LT.O.ODO)  AZ=AZ  -  180. ODO 
IFIN=60. ODO/DELTAU+1 . ODO 
TAU=TAUI 
C 

C  INITIALIZE  IMAGE  WINDOW  TO  MIN.  CALC  OF  SUNGLINT  AND  RESTRICT  CA 

C  TO  AREA  OF  INTEREST. 

cal  1  realgaeOh/platmin/,  platmin) 
call  realgaeOh /platmax/,  platmax) 
call  realgae(7h/plone/, plone) 
call  realgae(7h/plonw/, plonw) 
if (plone. It. plonw)  plone=plone+360. 0 
C 

DO  13  I  =  1 f IFIN 

C  ORGITAL  LONG  DIF  BETWEEN  SUN  AND  SATELLITE 
CSDOBG=DCOS( (SL0BLG-TAU)*PI/1B0. DO) 

CCC  write(6, ' (3(f8.2,2x) ) ' ) SNDTAU, CSDTAU, PHI 

CALL  SUBLTLN (SNDTAU, CSDTAU, PHI , ALN, AZ , PHI2, ALN2, AZ2) 

ETA =-40. ODO 

■JFIN=80.  ODO/DELETA+l .  ODO 
C 

C  INITIALIZE  CALL  ELOC 

JE  =  1 
C 

DO  15  J=1 , JFIN 

GETA=-ETA 

ETARD=ETA*PI/180. 

SNTHET = ( R+HT ) *DS I N ( ETARD ) / R 
THET  RD  =DAS I N ( SNTHET ) 

C  SATELLITE  ZENITH  ANGLE  AT  VP 

THET A  =DABS ( THET  RD* 1 80 . DO / P I ) 

C  GEOCENTRIC  ANGLE  BETWEEN  SATELLITE  AND  VP 
S I G  =THETRD -ET ARD 
SNS I G  =DS I N ( S I G ) 

CSS I G  =DCOS ( S I G ) 

CCC  WRITE  (6,6006)  SNSLOT, CSSLOT, CSDOBG, SNSIG, CSSIG, ZETA 

CALL  GEOMET ( SNSLOT , CSSLOT , CSDOBG , SNS I G , CSS I G , ZETA , ALPHA ) 

CCC  WRITE  (6,6006)  ALPHA, THETA, WIND,  RI 

C 

C  FIND  LAT  AND  LONG  FOR  POINT  BEING  CONSIDERED 

C 

CALL  ELOC ( JE , REAL ( ETA ) , REAL (HT) , REAL (PHI), REAL ( PHI 2 ) , 

+  REAL (ALN), REAL (ALN2), PLAT, PLON, IERR)  ! EARTH  LOCA 

C 

if (plon. It. plonw)  then 
plont«plon+360.0 
else 

piont  *plon 
end  if 

IF (PLATMAX. GE. PLAT  .AND.  PLATMIN. LE. PLAT)  THEN 
IF( PLONW. LE. PLONT  .AND.  PLONE. GE. PLONT)  THEN 
CALL  INTERP(W, PLAT, PLONT, PLATMAX, PLONW, WIND) 

CALL  OCGLNT(ZETA, THETA, ALPHA, DFLOAT (WIND) ,RI) 
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RA=100. 0D0*RI 
C 
C 

CCC  WRITE  (6,6006)  GETA, TAU, RI , TAUI , ORBINC, SNINC 

IF(RI.GE. (.00100) )  THEN 

QXX383XZ22X222223XX3X2X222222223222222222222S3232222=22222Z222X2222X22  =C 

CALL  PLACE( REAL (GETA) , REAL (TAU) , REAL (RI ) , 1 , REAL (TAUI ) , 

+  REAL (ORBINC) , REAL (SNINC) , REAL (WIND) , PLAT, PLON  ) 

ENDIF 
END  IF 
ENDIF 

C  ask  jim  about  the  placement  of  je=je+l  should  it  be  just  befor 

c  call  place  inside  if  stmts  or  here??? 

JE= JE+1 

15  ETA=ETA+DELETA 

PHI =PHI2 
AZ  =  AZ2 
ALN  =  ALN2 
13  TAU=TAU+DELTAU 

C 

CCC  WRITE  (6,6006)  GETA, TAU, RI , TAUI , ORBINC, SNINC 

C 

CALL  PLACE (REAL (GETA) , REAL (TAU) , REAL (RI) ,2, REAL (TAUI )  , 

+  REAL (ORBINC) , REAL (SNINC) , REAL (WIND)  ) 

11  CONTINUE 

C 

c 

c 

5000  FORMAT (2i3) 

5010  FORMAT (61 (13F6.  1/) ) 

C 

6000  FORMAT ( /  /  ’  DATA  =  ’,3A6/) 

6001  FORMAT ( '  SATELLITE  HEIGHT  *' , F6. 0, » ( KM) ,  SATELLITE' 

1»  INCLINATION  ANGLE  =  ',F6.1) 

6002  FORMAT ( *  SOLAR  POSITION:  DECLINATION  =',F6.1,',  SOLAR  ', 

1' LONGITUDE  =',F6.1,A1) 

6003  FORMAT ( *  EQUATOR  CROSSING:  TIME  =  »,A4,'Z,  LONGITUDE  *’ , 

1F6. 1, Al) 

6004  FORMAT ( ’  SOLAR  LONGITUDE  RELATIVE  TO  LONGITUDE  OF  EQUATOR  ', 

1» CROSSING  , F6. 1 , ’  EASTWARD') 

6005  FORMAT ( r  CHART  NO. » , 12, * ,  WIND  SPEED  , F6. 1 , ' (M/S) ' ) 

6006  FORMAT ('  CHECK  ',6F15.6) 

C 

7001  FORMAT (12) 

7002  FORMAT ( 12, '032100'  ) 

7003  FORMAT (A6,A4) 

7004  FORMAT (2F2.0) 

C 

9939  END 
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SUBROUT I NE  OCGL NT ( ZET A , THETA , GAL  PHA , W , R I ) 

IMPLICIT  REAL  *8(A-H),  REAL  *B(0-Z) 

SAVE 

C  CALCULATES  ANGULAR  REFLECTION  FROM  OCEAN  SURFACE 
C  INPUT  (ALL  IN  DEGREES) -ZETA“SOLAR  ZENITH  ANGLE, 

C  THETA=ZENITH  ANGLE  OF  REFLECTED  RAY,  GAL PHA = 

C  AZIMUTH  ANGLE  OF  REFLECTED  RAY  RELATIVE  TO  SUN(I.E., 

C  AZIMUTH  OF  SUN=0) .  GALPHA  MAY  BE  EXPRESSED  AS 

C  0-360, -180-+180,  ETC.  **IF  THETA =0,  GALPHA  MUST  BE  ZERO.** 
C  W«WIND  SPEED  IN  M/SEC 

C  OUTPUT-  RI “REFLECTED  INTENSITY  (/STERADIAN) , 

C  RELATIVE  TO  INCOMING  SOLAR  FLUX  OF  1 
C 

PI  =  DACOS  < - 1 . ODO ) 

HPI  =  PI/2. ODO 
RPD  *  PI/ 180. ODO 

ANGINC “ANGULAR  INCREMENT, IN  RADIANS,  USED  IN  CALCULATING 
SLOPE  PROBABILITY  AND  SOLID  ANGLE  OF  REFLECTED  BEAM 

ANGINC  =  0.01 DO 

STPNMB=NUMBER  OF  STEPS  PER  SIDE  USED  IN  INTEGRATION 
FOR  SOLID  ANGLE 

STPNMB  =  20. ODO 

IF  SUN  IS  BELOW  HORIZON,  SKIP  COMPUTATION 
IFCZETA.LT. 90. OGO  TO  8 
RI  =  O.ODO 
GO  TO  5 

8  SNZET  =  DSIN(ZETA*RPD) 

CSZET  =  DCOS ( ZETA*RPD) 

SNSQZT  =SNZET*SNZET 
SNTHET  =  DS I N ( THET A* RPD ) 

CSTHET  =  DACOS (THETA*RPD) 

SNSQTH “SNTHET *SNTHET 
C 

C  NEXT  2  STEPS  CONVERT  GALPHA  TO  ALPHA,  A  POSITIVE  ANGLE 
C  BETWEEN  0  AND  180  DEGREES.  (PERMISSABLE  SINCE  REFLECTION 
C  PATTERN  IS  SYMNETRICAL  ABOUT  0-180  LINE). 

AL  PHA “DABS  <  GA  L  PHA ) 

C 

IFCALPHA.GT. 180.0)  ALPHA  *  360. ODO -ALPHA 
CSALF  *  DCOS (ALPHA  *RPD) 

C  S I GSQ “VARIANCE  OF  THE  CIRCULAR  GAUSSIAN  DISTRIBUTION  OF 
C  THE  NORMALS  OF  THE  SLOPES  OF  THE  OCEAN  SURFACE 
SIGSQ  *  O . 005 1 2D0* W+0 . 003D0 
CS2BET “CSZET  *CSTHET  +SNZET *SNTHET *CSALF 
C 

C  CSBET “COSINE  OF  ANGLE  OF  INCIDENCE  (AND  ANGLE  OF  REFLECTION) 

C  NECESSARY  TO  REFLECT  INCOMING  RADIATION  IN  A  DIRECTION 
C  THETA, ALPHA 

CSBET  »  DSQRT (0. 5D0* ( 1 . 0D0+CS2BET) ) 

C  IF  CSBET  IS  NEGATIVE, SLOPE  IS  SHADOWED 


no  nonnnnnnn  n  n  nnnnnn  non  no 


IF(CSBET. GT. 0. 0)  GO  TO  9 
RI  =  O.ODO 
GO  TO  5 

RHO=REFLECTIQN  COEFFICIENT  (APPROXIMATION  TO  FRESNEL  LAW) 

9  RHO  =DEX  P ( -6 . 09D0*CSBET ) +0 . 0 1 9D0*CSBET 

SRT  =  DSQRT (SNSQZT+SNSQTH+2. 0*SNZET  *SNTHET*CSALF) 

TNTHN  IS  SEA  SURFACE  SLOPE  (TANGENT  OF  ZENITH  ANGLE  OF  SURFACE 
NORMAL)  REQUIRED  TO  REFLECT  AT  GIVEN  THETA  AND  ALPHA 
TNTHN  =  SRT / (CSZET  +  CSTHET) 

THTNC  *  DAT AN ( TNTHN) 

IF  THE  REFLECTED  DIRECTION  IS  THE  SPECULAR  POINT,  THE  SURFACE 
SLOPE  IS  ZERO,  AND  THE  SURFACE  NORMAL  IS  VERTICAL.  THE 
AZIMUTH  OF  THE  NORMAL,  ALFNC,  IS  UNDEFINED,  AND  SRT*0. 

ALFNC  IS  DEFINED  TO  BE  ZERO  IN  THIS  CASE  BY  STATEMENT  10, 
SINCE  ALPHA  FOR  THE  SPECULAR  POINT  IS  180  OR,  FOR  AN  OVERHEAD 
SUN, 0. 

IF (SRT. EQ. 0.0)  GO  TO  10 
CSALNC  =  (SNZET  +  SNTHET*CSALF ) /SRT 

10  IF(ALPHA. EQ. 0.0. OR. ALPHA. EQ. 180.0)  CSALNC  =  l.ODO 
IF(ALPHA.EQ. 180.0. OR. THETA. GT.ZETA)  CSALNC  =  -l.ODO 
IF (CSALNC. GT. 1.0)  CSALNC  =  l.ODO 
IF(CSALNC. LT. -1.0)  CSALNC  =  -l.ODO 

ALFNC  =  DACOS (CSALNC) 

NEXT  SECTION  SETS  UP  THE  ANGULAR  REGION  OVER  WHICH 
DETERMINATION  OF  SLOPE  PROBABILITY  AND  SOLID  ANGLE 
ARE  MADE 

L  AND  U  AT  END  OF  NAMES  REFER  TO  LOWER  AND  UPPER  LIMITS 
RESPECTIVELY.  THTN=ZENITH  ANGLE  OF  SLOPE  NORMAL,  AND 
ALFN  IS  ITS  AZIMUTH.  THTNL  AND  ALFNL  CANNOT  BE  NEGATIVE, 
THTNU  CANNOT  EXCEED  PI/2,  AND  ALFNU  CANNOT  EXCEED  PI. 

THTNL  =  THTNC  -  ANGINC 
IF (THTNL . LT. 0. 0)  THTNL  =  O.ODO 
THTNU  =  THTNC  +  ANGINC 
IF (THTNU. GT. HPI )  THTNU  =  HPI 
DTHTN  =  0.5D0-M  THTNU -THTNL) 

ALFNL  =  ALFNC  -  ANGINC 
IF ( ALFNL . LT. 0. 0)  ALFNL  *  O.ODO 
ALFNU  =  ALFNC  +  ANGINC 
IF (ALFNU. GT. PI)  ALFNU  =  PI 
DALFN  =  0.5  * (ALFNU -ALFNL ) 

TL  »  DT AN (THTNL) 

TU  =  DTAN( THTNU) 

P=PROBABIL ITY  OF  HAVING  SLOPES  WITHIN  THE  DEFINED  LIMITS 

P  =  (DEXP( -TL*TL/SIGSQ) -DEXP ( -TU*TU/SIGSQ) )  *DALFN/PI 

C  THE  REMAINDER  OF  THE  PROGRAM,  THROUGH  STATEMENT  3,  CALCULATES 
C  THE  SOLID  ANGLE  INTO  WHCIH  RADIATION  IS  REFLECTED  BY 
C  SLOPES  HAVING  NORMALS  LYING  BETWEEN  THTNL  AND  THTNU, 
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C  AND  ALFNL  AND  ALFNU.  **  THIS  SOLID  ANGLE  MUST  BE 
C  COMPUTED  TO  AGREE  AS  EXACTLY  AS  POSSIBLE  WITH  THE 
C  RANGE  OF  NORMALS,  SINCE  LARGE  ERRORS  IN  THE 
C  FINAL  RESULT  MAY  OTHERWISE  OCCUR.** 

C 

C  NSTP=NUMBER  OF  VALUES  PER  SIDE  FOR  WHICH  CALCULATED  (INCLUDES  POINT 

c  zero:> 

NSTP  =  STPNMB  +  l.ODO 

SPTHN  AND  SPALN  ARE  STEP  INTERVALS  OF  NORMAL  THETA  AND  ALPHA 
SPTHN  =  2. ODO*DTHTN/STPNMB 
SPALN  =  2. ODO*DALFN/STPNMB 
1421  FORMAT ('  HERE') 

SGN  =  l.ODO 

INTEGRATION  ALWAYS  STARTS  AT  SMALLEST  VALUES  OF  NORMAL  THETA  AND 
ALPHA,  INCREASES  ALPHA,  THEN  THETA,  DECREASES  ALPHA,  THEN 
THETA.  IN  THIS  WAY,  THE  PROPER  SIGNS  ARE  PRESERVED. 

SGN  CONTROLS  THE  INCREASE  OR  DECREASE  OF  THE  STEP  INTERVALS. 

THTST  =  THTNL 

ALFST  =  ALFNL  -  SPALN 

OMEGA =SOL ID  ANGLE  INTO  WHICH  RADIATION  IS  REFLECTED 
OMEGA  =  O.ODO 
CSTHST  =  DCOS( THTST) 

SNTHST  =  DSIN (THTST) 

DO  3  1=1,2 

SOME  CONDITIONS  OF  INTEGRATION  RESULT  IN  ZERO  CONTRIBUTION  TO 
THE  SOLID  ANGLE.  WHILE  THE  PRECISION  OF  THE  CALCULATIONS 
(DOUBLE  PRECISION  ON  IBM  360)  IS  SUFFICIENT  TO  MAINTAIN 
PRACTICAL  SIGNIFICANCE  IN  THE  SOLID  ANGLE,  IT  IS  POSSIBLE 
IN  CERTAIN  CASES  TO  COMPUTE  SINES  OR  COSINES  WHICH  ARE  SLIGHTLY 
GREATER  THAN  1.  SINCE  THESE  WILL  CAUSE  PROGRAM  FAILURE,  THE  CASES 
FOR  WHICH  THIS  MAY  OCCUR  ARE  TREATED  SEPARATELY,  AND  THE 
INTEGRATION  STEPS  ARE  OMITTED.  THESE  ARE  CONTROLLED  BY  THE  TWO 
GROUPS  OF  *IF*  STATEMENTS  BELOW,  THE  FIRST  GROUP  WHEN 
THE  NORMAL  ZENITH  ANGLE=0,  AND  THE  SECOND  GROUP 
WHEN  THE  NORMAL  AZITUTH  IS  0  OR  PI. 

IF (THTNL . EQ. 0. 0. AND. I . EQ. 1 )  ALFST  =  ALFNU 
IF (THTNL . EQ. 0. 0. AND. I . EQ. 1 )  CSALST  =DCOS(ALFST) 

IF (THTNL. EQ. 0.0. AND. I.EQ.  1)  THTST  =  THTNL 
IF (THTNL . EQ. 0. 0. AND. I . EQ. 1 )  GO  TO  11 

DO  1  J=1 , NSTP 

ALFST  =  ALFST  +  SGN*SPALN 

CSALST  =  DCOS( ALFST) 

IF(I.EQ. 1. AND. J.EQ. 1. AND. ALFNL. EQ. 0.0)  CSALST  =  l.ODO 
IF(I. EQ. 2. AND. J.EQ. NSTP. AND. ALFNL. EQ. 0.0)  CSALST  =  l.ODO 
IF ( I. EQ. 1. AND. J.EQ. NSTP. AND. ALFNU. EQ. PI )  CSALST  =  -l.ODO 
IF( I . EQ. 2. AND. J.EQ. 1 . AND. ALFNU. EQ. PI )  CSALST  =  -l.ODO 

SLDANG  CALCULATES  THE  INCREMENT  OF  SOLID  ANGLE 
ccc  PRINT  *, ’  CHECK  SLDAN61  ' , snzet , csz*t , 


cn 


ccc  + 


xnthst, csthst, csalst, j, omega 
CALL  SLDANG (SNZET, CSZET , SNTHST , "CSTHST , CSALST f  J , OMEGA > 
ccc  PRINT  CHECK  SLDANG1  ' , snzet , cszet , 

ccc  +  xnthst, csthst, csalst, j, omega 

1  CONTINUE 
C 

IF  ( ALFNL .  EQ.  0.  0.  AND.  I .  EO.  2)  GO  TO  4 
IF  CALFNU. EQ. PI . AND. I . EQ. 1 )  THTST  =  THTNU 
IF (ALFNU. EQ. PI . AND. I . EQ. 1 )  SNTHST  =  DSINC THTST) 
IF<ALFNU. EQ. PI . AND. I . EQ. 1 )  CSTHST  =  DCOS(THTST) 

IF  CALFNU. EQ. PI . AND. I . EQ. 1 )  GO  TO  21 
11  THTST  =  THTST-  SGN*SPTHN 
DO  2  J=1 , NSTP 
THTST  =  THTST  +  SGN*SPTHN 
SNTHST  =  DS IN (THTST) 

CSTHST  =  DCOSC THTST) 

ccc  PRINT  *, f  CHECK  SLDANG2  ’, snzet , cszet , 
ccc  +  xnthst, csthst, csalst, j, omega 

CALL  SLDANG (SNZET, CSZET , SNTHST , CSTHST , CSALST , J, OMEGA) 
ccc  PRINT  CHECK  SLDANG1  snzet , cszet , 

ccc  +  xnthst, csthst, csalst, j, omega 

2  CONTINUE 

21  ALFST  =  ALFNU  +  SPALN 

3  SGN  =  -1.0 
RI  =  RHO*P*CSBET /OMEGA 
RETURN 
END 
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SUBROUT I NE  SL DANG  ( SNZET , CSZET , SNTHST f  CSTHST , CSAL  ST , J ,  OMEGA ) 
C  SLDANG  CALCULATED  THE  INCRENENT  OF  SOLID  ANGLE 
C  OF  REFLECTED  ENERGY,  OMEGA.  SEE  OCGLNT  FOR  DISCUSSION 

C  OF  DEFINITIONS.  THE  INTEGRATION  IS  PERFORMED  ON  A 

C  LAMBERT  AZIMUTHAL  EQUAL  AREA  PROJECTION  OF  THE 
C  UPWARD  HEMISPHERE,  CENTERED  ON  THE  SOLAR  DIRECTION 
C  (I.E.  ,  THETA=ZETA,  ALPHA*0) . 

C 

C  2*BETA=THE  ANGLE  BETWEEN  THE  SUN  AND  THE  REFLECTED  RAY. 
IMPLICIT  REAL  *B(A-H),  REAL  *8CO-Z> 

SAVE  STRB, STRM 

CSBT  =  CSZET  *CSTHST  +  SNZET*SNTHST*CSALST 
C  2*SIN(BETA)  =2*SNBT “RADIUS  OF  REFLECTED  RAY  ON  PROJECTION 
SNBT  =  DSQRT  (  1 . 0-CSBT*CSBT) 

C 

I F  ( SNBT . NE .0.0)  GO  TO  2 
CSMU  =  CSAL ST 
GO  TO  3 

AMU=AZ IMUTH  AT  PROJECTION  CENTER  FROM  DIRECTION  TO  VERTICAL 
(IF  SUN  IS  VERTICAL,  AMU  MEASURED  FROM  LINE  ALPHA=18Q) 

2  CSMU  =  (SNZET *CSTHST  -  CSZET*SNTHST*CSALST) /SNBT 
*IF*  STATEMENTS  FORCE  CSMU  TO  BE  1  OR-1  IN  APPROPRIATE 
CASES,  IN  ORDER  TO  AVOID  PROGRAM  ERROR  DUE  TO  ROUNDOFF 
WHEN  ARC  COS  IS  TAKEN 

3  IFCCSALST.EQ. -1.0)  CSMU  =  l.ODO 

IF(CSALST.EQ. 1.0. AND. SNTHST. LE. SNZET)  CSMU  =  l.ODO 
IFCCSALST.EQ. 1.0. AND. SNTHST. GT. SNZET)  CSMU  =  -l.ODO 
AMU  =  DACOS (CSMU) 

IF  J-l,  THE  FIRST  POINT  ON  A  SIDE  HAS  BEEN  COMPUTED,  AND 
AN  INCREMENT  OF  AREA  CANNOT  BE  FOUND.  RETURN 
AND  COMPUTE  FOR  2ND  POINT 

IFCJ.EQ.l)  GO  TO  1 

STRB=SNBT  FOR  J-l,  STRM=AMU  FOR  J-l 

IF  1ST  POINT  COMPUTED  FOR  SIDE  WAS  AT  ORIGIN  (OCCURS  FOR  SUN 
AT  OR  NEAR  ZENITH)  FORCE  AZIMUTH  DIFFERENCE  TO  BE  ZERO 
BY  SETTING  AMU( 1 ) =AMU(2) .  THE  GENERAL  EQUATION  MAY 
RESULT  IN  A  SUBSTANTIAL  ERROR  FOR  AMU(l),  SINCE  IT  IS 
REALLY  UNDEFINED 

IF (STRB. EQ. 0.0)  STRM  *  AMU 
RAD “MEAN  RADIUS  FOR  INCREMENT 
RAD  =  SNBT  +  STRB 

DMU-AZ IMUTH  DIFFERENCE  FOR  INCREMENT 
DMU  “  AMU  -  STRM 

ADD  INCREMENT  OF  SOLID  ANGLE  TO  PREVIOUS  SUM 
OMEGA  *  OMEGA  +0. SDO*RAD*RAD*DMU 
STORE  SNBT  AND  AMU  FOR  USE  IN  NEXT  INCREMENT 

1  STRB  *  SNBT 

STRM  =  AMU 
RETURN 

END  A-10 


SUBROUTINE  GEOMET (SNSLOT, CSSLOT, CSDOBG, SNSIG, CSSIG. ZETA  ALPHA 
IMPLICIT  REAL  *8(A-H),  REAL  *8(0-Z)  f 

C  PREFI XED-SN=SINE,  CS=COSINE 

C  SIG=GEOCENTRIC  ANGLE  OF  VIEWED  POINT  (VP)  FROM  ORBITAL  TRACE 
C  ZETA=SOLAR  ZENITH  ANGLE  AT  VP 

C  ALPHA*AZ IMUTH  AT  VP  FROM  SUN  TO  SATELLITE  (GE. 0, LE. 180) 

PI =DACOS< -l.ODO) 

RTD=180. ODO/PI 

CS  ZET “SNSLOT *SNS I G+CSSLOT*CSS I G*CSDOBG 
SNZET  =DSQRT ( 1 . ODO -CS  Z  ET*CSZET ) 

ZETA  =DACOS  <  CSZET ) *RTD 

CSALF  =  (SNSLOT  *CSSIG-CSSLOT*SNSIG*CSDOBG) /SNZET 
IF (SNSIG. NE. 0. 0) GO  TO  4 
CSALF= l.ODO 
GO  TO  5 

4  CSALF=-CSALF*SNSIG/ABS(SNSIG) 

5  if (dabs( csalf ) . gt .  l.OdO)  then 

print  greater  than  one...' 
csalf  =  dint(csalf) 
end  if 

AL  PHA  =DACOS ( CSALF ) *RTD 

RETURN 

END 


SUBROUTINE  SLORBCCSNINC, CSINC, SNDEC, CSDEC, SNLMDA, CSLMDA, SNSLOT, 
1CSSLOT , SLOBLG) 

IMPLICIT  REAL  *8<A-H>,  REAL  *8<0-Z> 

SAVE 

C  PREFIa'ES-SN*SINE,CS=COSINE 

C  INC “ORBITAL  INCLINATION  (LE.90),  DEC=SOLAR  DECLINATION, 

C  LMDA=SOLAR  LONGITUDE  REFERRED  TO  ASCENDING  NODE, 

C  <0-360  DEG, W) ,  SLOT=SOLAR  ORBITAL  LATITUDE,  SLOBLG=SOL AR 
C  ORBITAL  LONGITUDE. 

C 

RTD  =  180. O/DACOS  < - 1 . ODO  > 

SNSLOT =CS I NC*SNDEC -SN I NC*CSDEC*SNLMDA 
CSSLOT  =DSQRT < 1 . ODO -SNSLOT*SNSLOT> 

CSSLOG=CSDEC*CSLMDA/CSSLOT 

SNSLOG=  <SNINC*SNDEC+CSINC*CSDEC*SNLMDA) /CSSLOT 

SLOBLG =DACOS<CSSLOG)*RTD 

IF (SNSLOG. LT. 0. 0) SL0BLG«360. ODO -SLOBLG 

RETURN 

END 
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SUBROUTINE  SUBLTLN(SNDTAU,  CSDTAU,  PHI ,  ALN,  AZ,  PHI2,  ALN2,  AZ2) 

IMPLICIT  REAL  *8(A-H),  REAL  *S(0-Z) 

PI »DACOS< -1 . ODO) 

D2R  =  PI/180.0D0 
R2D  *  180. ODO/PI 
PHIR  =  PHI*D2R 
AZR  *  AZ*D2R 
CSPHI  =  DCOS(PHIR) 

SNPHI  *  DSIN(PHIR) 

PH 1 2  =  DAS INC  CSDT AU*SN  PH I  +  SNDTAU*CSPHI*DCOS ( AZR) ) 

DELLN  =  DACOSC (CSDTAU  -  SNPHI*DSIN  ( PHI2)  )  /  (CSPHI*DCOS  (PHI2) )  )  *R2D 
ccc  AZ2  *  DASIN(CSPHI*DSIN(AZR)/DC0S(PHI2))*R2D+quad 
AZ2  =  CSPHI*DSIN(AZR)/DC0S(PHI2) 

C***********CHECK  FOR  QUADRANT***************************************** 
if(  (az.ge.O.OdO  .and.  az . le. 90. OdO)  .or. 

*  (az. gt. 270. OdO  .and.  az. It. 360. OdO)  )  then 

az2sdasin(az2)*r2d 
else 

az2  =  180. OdO  -  dasin Caz2)*r2d 
end  if 

C*********************************************************************** 
PH I 2  =  PHI2*R2D 
if (phi. It. phi2)  then 
ALN2  *  ALN  -  DELLN 
else 

ALN2  *  ALN  +  DELLN 
endif 
RETURN 
END 


subroutine  initgae 


c  written  by  a. b.  caughey  18  aug  1988 

c 

c  subroutine  to  access  gae  parameter  block 

c 

c 

integer  xnu 11 string, xprmdim, xcont, xsuccess, xci,xcs,xsave 
c  include  'subdef.fin' 

parameter  (xprmdim=1000, xsuccess=l, xcont=2) 
integer  block(xprmdim) 
c 

character  cmsg*72,  creply*72 
integer  imsg(18) 
integer  reply<18) 

equivalence  C  cmsg,  imsg) ,  C  creply,  reply) 
c 

parameter  (max_word =1024) 
integer*2  ibuf (max_word) 
integer  cmdblk(80) 

integer  factor, r len, band, type, x, name (6) , strg(7) 
character  cname*24, cstrg*28 
equivalence  (cname,name) 
real  y 

integer  idsrn,odsrn 
integer  output 
logical  ok 
c 

integer  cnm(3) , unitno 
character*28  blknam 
integer  ipackos(.'8) 
save  block 


c 

c 

c 


c 

c 

c 


***get  parameter  block  from  tm 

call  xrinim( block, xprmdim, xcont, istatus) 
if (istatus. ne. xsuc cess)  then 
13  format (' xrinim  error,  status  =  ',i6) 
wr ite< cmsg, 13)  istatus 
call  printpC 1 , imsg, 72, 0,0,0, 0, 0,0, 0,0) 
goto  1000 
endif 
return 

•**get  parameter  block  from  disk  file 


entry  initblk (blknam) 
lun=12 


call  xsch2p(bl knam,  ipackos) 

call  xr rdb ( ipackos, lun, block, xprmdim,  1 ,  istatus 
if (istatus. ne. xsuccess)  then 

write(cmsg, ' ("xrrdb  error,  status  ■  ",i6)’) 
call  printp( 1 , imsg, 72, 0,0, 0,0,0, 0,0,0) 
goto  1000 

endif  A-14 


) 

istatus 


return 


c  ***write  parameter  block  to  disk  file 
c 

entry  wrtblk (blknam) 

1 un=12 

call  xsch2p(blknam, ipackos) 

call  xqwrtb( ipackos, lun, block, istatus) 

if (istatus. ne. xsuccess)  then 

write(cmsg, ’ ("xqwrtb  error,  status  =  ",i6)’)  istatus 
call  printpC 1 , imsg, 72, 0,0, 0,0, 0,0, 0,0) 
goto  1000 
end  if 
return 
c 

c  **•*  open  i/o  files 
c 

entry  f i leop ( cnm, unitno) 

call  xrstrp(block,cnm,10, name, istatus) 

if ( istatus. ne. xsuc cess )  then 

21  format (’ xrstrp  error:  parameter  *  name,  status  =  ’,i6) 
write ( cmsg, 21 )  istatus 

call  printpCl,  imsg, 72, 0, 0, 0, 0, 0, 0, 0, 0) 
goto  1000 
end  if 

inquireCf ile=cname, exist =ok) 

if (.not. ok  .and.  (unitno. eq. 5  .or.  unitno. eq. 1 ) )  then 

31  format (’ input  file  does  not  exist’) 
write( cmsg, 31) 

call  printp< 1 , imsg, 25, 0, 0,0, 0, 0, 0, 0, 0) 
goto  1000 
endif 

open(unit=unitno, f ile=cname) 

C 

return 

c 

c  #**  get  real  parameter 
c 

entry  realgae( cnm, xx) 

call  xrreal (block, cnm, l,xx, istatus) 

if (istatus. ne. xsuccess)  then 

13  format (’ xrreal  error:  parameter  *  x,  status  =  ’,i6) 
wr ite( cmsg, 19)  istatus 
call  printp( 1 , imsg, 72,0,0,0,0,0,0,0,0) 
goto  1000 
endif 
return 
c 

c  ***  get  integer  paramter 
c 

entry  intgae(cnm, ix) 

call  xrintg(block,cnm,l,ix, istatus) 

if (istatus. ne. xsuccess)  then 

17  format (’ xrintg  error:  parameter  =  ix,  status  =  ’,i6> 
write( cmsg, 17)  istatus 
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call  printpCl,  imsg,  72,0,0,0,0,0,0,0,0) 
goto  1000 
endif 
return 
c 

c  ***  get  string  paramter 
c 

entry  strgae ( cnm, cstrg) 

call  xrstrp(block,cnm,7,strg, i status) 

if (istatus. ne. xsuccess)  then 

writeCcmsg, * ("xrstrp  errors  parameter  =  M,3a4,”  status  =  ",iS)' 
+  CcnmC i) , i=l , 3) , istatus 

call  printpCl, imsg, 72, 0,0, 0,0, 0,0, 0,0) 
goto  1000 
endif 

call  xsp2ch (strg, cstrg) 
return 
c 

c  ***  write  a  real  parameter  to  parameter  block 
c 

entry  wrtreal < cnm, xx ) 

call  xqreal(block,cnm,l,xx,103, istatus) 
if (istatus. ne. xsuccess)  then 

writeCcmsg, ' ("xqreal  error:  parameter  =  xx,  status  =  ",i6)') 

+  istatus 

call  printpCl, imsg, 72,0,0, 0, 0,0, 0,0,0) 
goto  1000 
endif 
retur n 
c 

c  **■*  write  a  integer  parameter  to  parameter  block 
c 

entry  wr iteint ( cnm, ix ) 

call  xqintgCblock, cnm, 1, ix, 103, istatus) 
if ( istatus. ne. xsuc cess)  then 

write( cmsg, ’ < " xqintg  error:  parameter-  ",3a4,",  status=  ",i6)r) 
+  ( cnm< i) , i=l , 3) , istatus 

call  printpCl, imsg, 72, 0, 0, 0, 0, 0, 0, 0,  0) 
goto  1000 
endif 
return 
1000  stop 
end 
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CCC*optic>n  trace  on 

SUBROUTINE  PLACE ( X , Y, RI , KEY, TAUI , ORBINC, SNINC, W, PLAT, PLON) 

SAVE 

C 

COMMON  /ABC/  I DATE (3) , I PLOT, L (5) 

C 

DIMENSION  VAL (5, SOO, 4) , XB(4) , XTC4) 

CHARACTER*12  reflfile  !  output  file  for  reflectance 

DATA  LL , ICNT / 1 , 0/ 

C  MAXIMUM  SIZE  OF  OUTPUT  ARRAY 

DATA  MXMAX , NYMAX / 61, 17/ 

C 

C 

I F  <  L  L . EQ . 0 ) GO  TO  101 
C 

XT  ( 1 ) =TAUI 
XB<1) *TAUI+15.0 
DO  100  LP=2, 4 
LM=LP-1 

XT  CLP) =XT(LM)+15.0 

100  XBCLP) =XBCLM)+15.0 

c 

c 

LL  =0 

101  CONTINUE 

IF  CABS ( X ) .GT. 40.0)  RETURN 
IF  <  KEY . EQ. 2)  GO  TO  4 
ICNT  =ICNT+1 
LCMP=5 

IF(Y.GE.XTCl). AND. Y. LE. XB< 1 ) )  LCMP=1 

IF  < Y. GT . XT<2) . AND. Y.LE.XB<2)>  LCMP=2 

IF(Y.GT.XT(3). AND. Y. LE. XB<3) )  LCMP=3 

IFCY. GT. XT ( 4 > . AND. Y. LE.  XB<4) )  LCMP=4 

IF  (LCMP.  EG).  5)  RETURN 

IF(L (LCMP). GE. BOO) RETURN 

LCLCMP!)  =  L  (LCMP)  +  1 

VAL (1,L< LCMP), LCMP)  =  Y 

VAL (2, L (LCMP) , LCMP)  =  X+23 

VAL (3, L (LCMP), LCMP)  =  RI*100 

VAL  <  4,  L  (LCMP) ,  LCMP!)  =  PLAT*3600. 

VAL (5, L( LCMP) , LCMP)  =  PL0N*3600. 

RETURN 

C 

C 

4  XMaO. 

DO  16  KK*1 , 4 

IF < L ( KK) . EQ. 0) GO  TO  16 

XM*XM+1. 

KL  =  5  -  KK 
16  CONTINUE 

C 
C 

C  PLOT  THE  VALUES 
C 

MM=XM 
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nnn  wn  n  n  nn  ono  no 


C 


SI ZE=0. 14 

output  file  opened  in  subroutine  initgae 
call  f i leop( lOh/ref If i le/ , 9) 
rewind  9 
I COUNT =1 


DO  3£  KK  =  1 ,  MM 
NUMV=L ( KK) 

XK=KK-1 

A=XK*4. 546875+1.0 
DO  36  1=1, NUMV 

CONVERTS  65.680  TO  0.  VALUES  IN  X  DIRECTION  TO  1.  TO  4. 546875* XK+1 . 
POSITION  IN  X  DIRECTION  =  ( VAL  Cl ,  I ,  KK)  -XT  <  KK) )  *4. 3786*0. 06922769+A 
XNUM= ( VAL  ( 1 , I , KK)  -XT(KK) )*0. 30312+A 

CONVERTS  -40  TO  +40  VALUES  IN  Y  DIRECTION 
YNUM= ( VAL (2,1, KK) +40. ) *. 096875  +  1. 

FVAL  =VAL  <3, I , KK) *10. 

MX  =  <  XNUM-1 . ) /0. 303125  +  1.2 
NY=(YNUM-1. )/0. 484375  +  1.2 
LAT  =NINT (VAL (4, I , KK) ) 

LON =N I NT (VAL  <5, I , KK) )  -  39600 

IF (MX . GE. 1  .OR.  MX. LE.MXMAX  .OR.  NY. GE. 1  .OR.  NY. LE. NYMAX ) then 
print  *,  'we  are  finally  printing  somethinq  ’ 

WRITE  <9,  ’  ("point'  ' ,  13. 3, 2<1X,  13)  ,  2  ( IX,  17)  ^3  ( 1 X ,  16.) )  '  ) 

+  I COUNT, MX, NY, LAT, LON,  REAL (FVAL) ,  real (fval) , i count 

endif 

I COUNT  =  I COUNT + 1 
IF( IPLOT.LT. 1)  GO  TO  36 

I CHAR  =FVAL  +0 . 5 
FVAL  =  ICHAR 
ICHAR=1 

IF (FVAL . GE.  9.5)  I CHAR =2 
IF (FVAL. GE.  99.5)  ICHAR=3 
IF (FVAL . GE. 999. 5)  I CHAR =4 
XSHFT  =SI ZE*0. 5 
YSHFT  =SIZE*ICHAR#0.84*0. 5 

CONTINUE 
CLOSE  (9) 

CLEAR  OUTPUT  ARRAY 

DO  3000  IX =1,3 

DO  2000  IY=1 , 800 
DO  1000  IZ  =  1, 4 

VAL (IX, IY, IZ)=O.ODO 
1000  CONTINUE 

2000  CONTINUE 
3000  CONTINUE 
RETURN 

C 

6000  FORMAT (6H  ****#, 'NOT  ENOUGH  POINTS  TO  PLOT' , 5H*****) 

C 


END 
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ccc*option  trace  on 

SUBROUTINE  ELOC( IEPH, SCANG, DALT, PHI , PHI2, ALN, ALN2, 

+  PLAT, PLON, IERR)  ! EARTH  LOCATE 

SAVE 

C  (FORMERLY  NAMED  ELOCAT) 

C  ELOCAT  GIVEN  THE  SCAN  TIME,  THE  SCAN  ANGLE,  AND  THE  THREE  ATTITU 

C  ANGLES,  COMPUTES  THE  LATITUDE  AND  LONGITUDE  OF  THE  PRINCIPAL  PO 

C  THE  ALTITUDE  OF  THE  SATELLITE  ABOVE  THE  EARTH,  AND  THE  ZENITH  A 

C  ABOVE  THE  PRINCIPAL  POINT. 

C  ELOCAT  WAS  BASED  ON  ROUGHLY  THE  SECOND  HALF  OF  THE  PROGRAM  CELCA 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c= 

c 


c 

c= 

c 

c 


FROM  M.  CHALFANT  OF  NOAA/NESS 
MODIFICATIONS  BY  G.  COLLINS,  ODSI  APRIL  1975 
MODIFICATIONS  BY  STEVE  REYNOLDS,  TELOS  COMPUTING,  NOV.  1977 
MODIFIED  BY  DONALD  C.  SCHERTZ  1978-1982  NEPRF  MONTEREY, CAL  IF 
INPUT 

IEPH  -  INDEX  INTO  EPHEMERIS  TABLES  FOR  SUBTRACK  POINT 
ISP  -  INDEX  INTO  SCAN  POINT  TABLES  FOR  SCAN  ANGLE 
OUTPUT 

PLAT  -  LATITUDE  OF  PRINCIPAL  POINT  (PRINCIPAL  POINT  IS 

THE  ASSUMED  CENTER  OF  THE  SCAN  AREA  ON  THE  EARTH. 

IN  DEGREES,  +  NORTHERN  HEMISPHERE,  -  SOUTHERN) 

PLON  -  LONGITUDE  OF  PRINCIPAL  POINT.  (DEGREES,  +  EASTWARD 
FROM  GREENWICH) 

IERR  -  ERROR  FLAG 

=0  -  NO  ERRORS 

*2  -  NEW  EPHEMERIS  TABLE  FOR  NEW  DAY  NEEDED 
*3  -  SPOT  OFF  EARTH 
=4  -  BAD  ORBITAL  DATA 
=5  -  UNKNOWN  TYPE  ERROR 

COMMON  /EPHE2/  CONTAINS  THE  EPHEMERIS  CALCULATIONS 
NEPH  -  NUMBER  OF  EPHEMERIS  POINTS 

STIME  -  STARTING  TIME  OF  EPHEMERIS  DATA  (SECONDS  AFTER  EPOCH) 

FT I ME  -  ENDING  TIME  OF  EPHEMERIS  DATA  (SECONDS  AFTER  EPOCH) 

ETIM  -  ELAPSED  SECONDS  SINCE  EPOCH  TIME 
ELAT  -  SUBSATELLITE  LATITUDE,  DEGREES,  +  NORTHERN 
HEMISPHERE,  -  SOUTHERN. 

ELON  -  SUBSATELLITE  LONGITUDE,  DEGREES,  +  EASTWARD 
FROM  GREENWICH  RANGING  0  TO  360 
EALT  -  SATELLITE  ALTITUDE  ABOVE  ELLIPSOID  (KM) 

(THE  LAST  FOUR  VALUES  ARE  PROVIDED  BY 
CALLS  TO  DATA  PAGING  ROUTINES) 


PARAMETER  (NSP=17) 
P I  * ACOS ( - 1 . ) 
PI0V2-PI/2.0 
TWOPI *2.0  #PI 
RADDEG=180. O/PI 
DEGRAD “PI/180.0 


33ZZ2S2S333SS3323B33333SB3S333S8ZSSS38SZS33SSB323BS2S3S3333333SSS  sQ 

IERR-0  UNIT  RETURN  ERR  INDICATOR 

SET  SCANG  FROM  SCAN  ANGLE  TABLE 

CHECK  IF  SUBTRACK  POINT  IS  SAME  AS  LAST  CALL, 
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C  SKIP  CALCULATION  OF  TRACKLINE  AZIMUTH  (A)  IF  SO. 

IF ( IEPH. GT. 1  )GO  TO  80 

C - - - C 

40  CONTINUE 

SPLAT =PHI*DEGRAD  ! CONVERT  TO  RADIANS,  1  =  ELAT 

SPLON  =ALN*DEGRAD  !2  =  ELON 

SPALT aDALT  ! SATELLITE  ALTITUDE  IN  KM,  3  =  EALT 

C  CHECK  FOR  QUICK  BYPASS  IF  SCAN  IS  ALONG  NADIR  (I.E.,  SUM  OF  SCAN 

C  AND  ROLL  ANGLES  IS  ZERO  AND  PITCH  ANGLE  IS  ZERO. ) 

IF (SCANG. EQ. 0. 0) GO  TO  280 

C  PICK  UP  NEXT  POINT  FOR  TRACKLINE  AZIMUTH  CALCULATIONS 

SL0N2*ALN2*DEGRAD  !2  =  ELON 

SLAT2=PHI2*DEGRAD  !1  *  ELAT 

C  THE  NEXT  SECTION  SOLVES  THE  SPHERICAL  TRIANGLE 

C  BETWEEN  THE  SUBPOINT,  A  FUTURE  SUBPOINT,  AND  THE  NORTH  POLE 

C  TO  DETERMINE  THE  TRACKLINE  AZ IMUTH(A) . 

SA=PI0V2-SLAT2 
CSA=COS (SA) 

SSA=SIN (SA) 

ccc  print  *,  ’  piov2  =  ’ ,piov2,’  splat  =  ’, splat,’  sb  =  ',sb 

SB=PI0V2 -SPLAT 
AC =SL0N2 -SPLON 
CSB=COS (SB) 

ccc  print  *,  ’  sb  =  ',sb 

SSB=SIN (SB) 

CAC=COS(AC) 


CKON  =SSA*SSB*CAC 
CSAB=CSA*CSB 
AKON  =CSAB+CKON 
SC=ACOS(AKON) 

IF (SC. EQ. 0. ) GO  TO  350 
CSC=COS(SC) 

SSC  =S I N ( SC ) 

ACSB  =CSB*CSC 
ASBC  =SSB*SSC 

ccc  print  *,  ’  csa  =  ’,csa,’  acsb  =  ’,acsb,’  asbc  =  ’,asbc 

ccc  print  *,  '  ssb  =  ’,ssb,’  ssc  *  ',ssc 

A* (CSA -ACSB ) / ASBC 
IF ( A. GE. 1. )G0  TO  350 
IF  (A.  LE.  -1.  )G0  TO  350 

A  = ACOS ( A ) * RADDEG  !  *****THIS  WAS  ALTERED  PER  JIM  CLARK***** 

50  IF (A. LT. 360. ) GO  TO  60 

A 3 A -360. 

GO  TO  50 

60  IF  (A.  GE.  0.  )G0  TO  70 

A*A+360. 

GO  TO  60 

70  A*A*DEGRAD 

C 
C 
C 

80  SCANG  sSCANG*DEGRAD 

PANAD  =SCANG 

C  THIS  SECTION  CALCULATES  THE  AZIMUTH  FROM  THE  SUB-SATELLITE  POINT 

C  TO  THE  PRINCIPAL  POINT  (AA). 


I F  (  SCANG  )  90 , 280 ,  1 00 
90  T  ANAD  =  -PANAD 

AA=A-PI0V2 
GO  TO  110 

100  TANAD =PANAD 

AA=A+PI0V2 

110  IF(AA) 120, 140,  130 

120  AA=AA+TWOPI 

GO  TO  110 

130  IFCAA.  LT.TWOPDGO  TO  140 

AA=AA-TWOPI 
GO  TO  130 

C  CALCULATION  OF  THE  ZENITH  ANGLE. 

140  RK0N=.0001569612*SPALT+1. 

Z  ANG  = AS I N  <  RKON*S I N  C  T ANAD ) ) 

C  TEST  FOR  EARTH  TANGENCY 

CNAL  =ASIN  (l./RKON) 

ccc  print  *, '  cnal  =  ' , cnal, '  tanad  =  ’, tanad 

IFCCNAL . LE. TANAD) GO  TO  360 

C  THIS  SECTION  SOLVES  THE  SPHERICAL  TRIANGLE  BETWEEN  THE  SUBPOINT, 

C  THE  PRINCIPAL  POINT  AND  THE  NORTH  POLE.  REMEMBER  IF 

C  A  B+C=180  THEN  A= (SUPPLEMENT  OF  B>  -C. 

SIDEA=ZANG-TANAD 
IF (SPLAT. LE. 0. ) GO  TO  150 
CNOS  =0 . 

SIDEB=PI0V2-SPLAT 
GO  TO  160 
150  CN0S=1 . 

SIDEB=PI0V2+SPL AT 
160  IF  (SIDEB.LT. 0. )G0  TO  370 

IF (AA) 370, 320, 170 
170  IF(AA-PI) 180,290,200 

180  CAPC  =AA 

W0E=1 . 

I F ( CNOS ) 370 , 220 , 1 90 
190  CAPC=PI -CAPC 

GO  TO  220 
200  W0E=0. 

CAPC  =TWOP I -AA 
IF (CNOS) 370, 220, 210 
210  CAPC=AA-PI 

220  SINA=SIN(SIDEA) 

SINB=SIN (SIDES) 

COSC =COS (CAPC) 

TEMP2 =S I NB*S I NA*COSC 
COSB  =COS ( S I DEB ) 

COSA  =COS (SI DEA ) 

SI DEC  =ACOS ( TEMP2+C0SB*C0SA ) 

SINC=SIN(CAPC) 

SI NEC  =SIN(SIDEC) 

CAPA  * AS I N ( S I NA*S I NC /S I NEC ) 

C  SEE  IF  CAPA  IS  IN  THE  SECOND  QUADRANT. 

COST  =COSB*COS ( S I DEC ) 

I F ( COST . GE . COSA ) CAPA =P I -CAPA 

C  NOW  DETERMINE  PRINCIPAL  POINT  LATITUDE  AND  LONGITUDE. 
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PLAT  =PI0V2-SIDEC 
I F  <  CNOS ) 370 , 240  f  230 
230  PLAT = -PLAT 

240  IF (WOE) 370, 250,260 

250  PLON=SPLON+CAPA 

IF(PLON. GT. TWOPI ) PLON=PLON-TWOPI 
GO  TO  270 

260  PLON=SPLON-CAPA 

IF  <  PLON. LT.O. )PLON'PLON+TWOPI 
C  FINAL  CLEAN-UP  AND  CONVERSION 

270  PLAT  =PLAT*RADDEG 

PLON=PLON*RADDEG 
GO  TO  380  ! RETURN 
C  SPECIAL  CASES 

C  SPECIAL  CASE.  THE  TRUE  NADIR  ANGLE  IS  ZERO  AND  THE  PRINCIPAL 

C  POINT  AND  THE  SUBPOINT  COINCIDE. 

200  TANAD=0. 

ZANG  =0 . 

PLAT  =SPLAT 
PLON=SPLON 
GO  TO  270 

C  SPECIAL  CASE  (LINE  TRIANGLE). 

290  PLON=SPLON 

PLAT  =SPLAT-SIDEA+PI0V2 
IF ( PLAT . GE. 0. ) GO  TO  300 
TLAT  =PLAT+PI0V2 
IF(TLAT. LT.O. )TLAT  =  -TLAT 
PLAT  * -TLAT 
GO  TO  310 

300  PLAT  =PLAT-PI0V2 

GO  TO  270 

C  SPOT  IS  OVER  THE  POLE,  REVERSE  THE  LONGITUDE. 

310  PLON=SPLON+PI 

IF ( PLON. GT. TWOPI ) PLON=PLON-TWOPI 
GO  TO  270 

C  SPECIAL  CASE  (LINE  TRIANGLE). 

320  PLON  =SPLON 

PLAT  =SPLAT+SIDEA-PI0V2 
IF ( PLAT. GT. 0. )G0  TO  330 
PLAT  =PLAT+PI0V2 
GO  TO  270 

330  TLAT  =PLAT-PIQV2 

IF (TLAT. LT.O. ) TLAT = -TLAT 
PLAT  =TLAT 
GO  TO  310 

£3333333333333333333333  ERROR  RETURNS  «333S«3»3B3333S33333333S33S=S33 

C  EPHEMERIS  TABLES  FOR  NEW  DAY  NEEDED 

340  IERR=2 

C  RESET  INDEX  FOR  NEW  EPHEMERIS  TABLE  ON  NEXT  ENTRY 

GO  TO  380 

C  IMPOSSIBLE  SITUATION,  MUST  BE  CAUSED  BY  BAD  ORBITAL  DATA. 

350  I ERR =4 

GO  TO  380 

C  OFF  THE  EARTH  RETURN. 

360  IERR=3 
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c 

c 

370 

380 


SET  DUMMY  VALUES 
GO  TO  380 

MISCELLANEOUS  ERRORS 

IERR  =5 

RETURN 

END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  IDIFDTG  C 


C  C 
C  PURPOSE:  TO  COMPUTE  THE  TIME  DIFFERENCE  BETWEEN  TWO  DISPLAY  CODE  C 
C  LEFT-JUSTIFIED  DATE-TIME  GROUPS  (DTGS) .  C 
C  C 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  C 
C  INPUT  BY  THE  USER:  NONE  C 
C  C 

C  c 

C  OUTPUT:  IDIFDGT  C 
C  C 
C  THIS  SUBROUTINE  IS  CALLED  BY  PROGRAM  SUNGLINT  C 
C  C 
C  THIS  SUBROUTINE  CALLS:  NONE  C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  VARIABLES 
C  INPUT 
C  DTG1 

C 

C  DTG2 


TYPE  ALLOWED 

FIRST  LEFT -JUSTIFIED  DATE -TIME  GROUP  CHAR  OOOOOOOO  - 

99123123 

SECOND  LEFT-JUSTIFIED  DATE-TIME  GRP  CHAR  OOOOOOOO  - 

99123123 


C  OUTPUT  C 

C  IDIFDTG  TIME  DIFFERENCE  IN  HOURS  BETWEEN  DTG2  C 

C  AND  DTG1  INTEGER  0  -  876767  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  METHOD:  THIS  FUNCTION  CONVERTS  DTG  INTO  HOURS  SINCE  1900  C 

C  THEN  SUBSTRACTS  2ND  DTG  FROM  1ST  DTG.  THE  DIFFERENCE  C 

C  IS  IN  HOURS.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

FUNCTION  IDIFDTG <DTGlfDTG2) 

C 

CHARACTER  *12  CDTG(2),  DTG1 ,  DTG2 
INTEGER  IMOSC 12) f IHRSC2) 

C 

C  NUMBER  OF  DAYS  IN  EACH  MONTH 


DATA  IMOS/31, 28,  31, 30,  31, 30,  31, 31, 30, 31, 30,  31/ 

C 

C  CALCULATE  NUMBER  OF  HOURS  FROM  190phcNZ9}T64491uvLhL44hBQu 

X6zo ! jT3L444h83 IKk#QG2 
C 

DO  100  1*1,2 

READ(CDTGd)  ,9001)  IYR,  IMO,  IDAY,  I  HR 
IF (MOD (IYR,4).EQ.O)  THEN 


n  n  n 


C  ADD  IN  YEARS 

IHRS  <  I ) =1 YR*366*24 
I MOS (2) =29 
ELSE 

IHRSCI)=IYR*365#24 
IM0SC2) =28 
END  IF 

C  ADD  IN  MOS 

DO  50  J=1 , I MO 

IHRSCI) =IHRS<I)  +  I MOS ( J) *24 
50  CONTINUE 

C  ADD  IN  DAYS  AND  HOURS 

IHRSCI) =IHRS(I)  +  <IDAY*24)  +  I HR 
100  CONTINUE 


TIME  DIFFERENCE  IN  HOURS  BETWEEN  IDTG1  &  IDTG2 

IDIFDTG=IHRS(2) -IHRSCI) 

RETURN 
9001  FORMAT (412) 

END 


A-25 


SUBROUTINE  INTERPCUVFLD, WLAT, WLON, FLAT, FLON, GSPD) 


C 

C  UVFLD  —  U  FIELD  IN  UVFLD( I I , J  J, 1 >  M/S 

C  V  FIELD  IN  UVFLDC 1 1 ,  JJ, 2)  M/S 

C  WLAT  --  LAT  OF  POINT  TO  BE  INTERPOLATED  TO  DEGREES 

C  WLON  --  LONG  OF  POINT  TO  BE  INTERPOLATED  TO  DEGREES 

C  FLAT  --  MAXIMUM  LAT  OF  UV  GRID  FIELD  DEGREES 

C  FLON  --  MAXIMUN  LONG  OF  UV  GRID  FIELD  DEGREES 

C  GSPD  --  WIND  SPEED  M/S 

C 

REAL  UVFLD (73, 144,2), GOBS (2) 

PARAMETER  (NREP=2, GRDLONSP *2. 5, GRDLATSP *2. 5) 

PARAMETER  (DLT*2. 5, DLN»2. 5) 

C  DLT  &  GRDLATSP  And  DLN  &  GRDLONSP  should  be  equi valenced ! 1 1 1 

CCC  EQUIVALENCE  (DLT, GRDLATSP) , (DLN, GRDLONSP) 

C 

C  INTERPOLATE  TO  OBSERVATION  POINT 
C 

DX=MOD( WLAT, GRDLATSP) 

DY  =MOD ( WLON , GRDLONSP ) 

6LAT  =WLAT-DX+GRDLATSP 
GLON =WLON-DY+GRDLONSP 
A*ABS (DY/DLT) 

B=ABS(DX/DLN) 

I =ABS (GLAT -FLAT) /2. 5+1 
J  =ABS ( GLON -FLON ) /2 . 5+ 1 
DO  100  K-l ,  NREP 

VI J  *  UVFLD ( I , J,  K) 

VIJP  s  UVFLD(I, J+l , K) 

VIPJ  *  UVFLD ( 1+1 , J, K) 

VPP  a  UVFLD (1+1, J+l, K) 

C 

C  PERFORM  BILINEAR  INTERPOLATION 

C 

GOBS ( K )  =  VIJ  +CVIJP  -VI J) *A  +CVIPJ  -VIJ)*B 
«<  +(VPP  +VIJ  -VIJP  -VIPJ) *A*B 

100  CONTINUE 

GSPD  =  ( SORT ( GOBS ( 1 ) **2  +  G0BS(2)**2)) 

C 

RETURN 

END 


APPENDIX  B 


SUNGLINT  PROGRAM  RONNING  INSTRUCTIONS 


The  Sunglint  algorithm  is  implemented  in  Global  Imaging's 
executi ve, operating  on  NEPRF's  HP9020.  There  are  two  modes  in 
which  to  run  the  sunglint  algorithm.  The  first  mode  requires  an 
input  image  with  a  navigational  file.  The  second  mode  requires 
the  user  to  supply  the  required  navigation  parameters. 

To  plot  sunglint  isopleths  on  a  navigated  image,  from  Global 
Imaging's  executive,  type  Htutor  SUN".  The  following  parameters 
must  be  assigned  values: 


PARAMETERS 

IMAGE  Input  image  name. 

Type  *  string 
Length  *  12 
Count  =  1 
Default^  "NONE" 

Gfile  Flag  "V"  is  set  if  a  navigational 

file  is  available  for  the  input  image. 

Type  *  string 
Length  *  1 
Default*  "Y" 

Valid  *  C "Y“, "N" > 

WIND  Wind  field  file  name. 

Type  *  string 
Length  *  12 
Count  *  1 

Default*  " INPUT. DAT" 


To  calculate  reflectance  for  a  given  area  and  time,  type 
"tutor  SUN".  The  output  will  be  a  file  of  reflectance  values  and 
their  location.  Values  must  be  assigned  to  the  following 
parameters: 


PARAMETERS 

IMAGE  Input  image  name. 

Type  »  string 
Length  *  12 


Count  =  1 
Default®  "NONE" 

Qfile  Flag  "Y"  is  set  if  a  navigational 

file  is  available  for  the  input  image. 

Type  =  string 
Length  =  1 
Default®  "Y" 

Valid  ®  ( "Y" , "N" ) 

YEAR  Year  of  image. 

Type  =  integer 
Count  =  1 
Default®  0 

MONTH  Month  of  image. 

Type  =  integer 
Count  =  1 
Default®  0 

DAY  Day  of  image. 

Type  =  integer 
Count  =  1 
Default®  0 

HOUR  Hour  of  image. 

Type  ®  integer 
Count  =  1 
Default®  0 

MIN  Minute  of  image. 

Type  =  integer 
Count  =  1 
Default®  0 

EQLONG  Longitude  of  Equator  crossing. 

Type  =  real 
Count  =  1 
Default®  0 

EQHEM  Hemisphere  of  Equator  crossing. 

Type  ®  string 
Length  ®  1 
Count  ®  1 
Default®  "W" 

SATINCL  Satellite  inclination  angle. 
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SATALT 


WIND 


pi atmi n 


pi atmax 


pi  one 


plonw 


V 


Type  =  real 
Count  =  1 
Default*  98.7 

Average  satellite  height  in  km. 

Type  -  real 
Count  =  1 
Default*  0 
Count*  1 

Wind  field  file  name. 

Type  =  string 
Length  =  12 
Count  =  1 

Default*  "INPUT. DAT" 


Latitude  of  southern  most  point  on  image. 

Type  =  real 
Default  *  12.5 
Count  =  1 

Latitude  of  northern  most  point  on  image. 

Type  =  real 
Default  *  30.0 
Count  *  1 

Longitude  of  eastern  most  point  on  image. 

Type  =  real 
Default  =  227.5 
Count  =  1 

Longitude  of  western  most  point  on  image. 

Type  *  real 
Default  =  252.5 
Coun  t  *  1 
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